Mon. Not. R. Astron. Soc. OOO.[THT2lf0000') Printed 12th October 2011 (MN WI^ style file v2.2) 

Analytic solutions to the accretion of a rotating finite cloud 
towards a central object - II. Schwarzschild spacetime 

Emilio Tejeda^*, Sergio Mendoza^* and John C. Miller^'^* 

1 SISSA & INFN, Via Bonomea 265, 34136, Trieste, Italy 

^ Instituto de Astronomia, Universidad Nacional Autdnoma de Mexico, AP 70-264, Distrito Federal 04510, Mexico 

^ Department of Physics (Astrophysics), University of Oxford, Keble Road, Oxford OXl 3RH, UK 



CN ■ 12th October 2011 



o 
O 



o 



% 



ABSTRACT 

Wc construct a general relativistic model for the accretion flow of a rotating finite cloud 
of non-interacting particles infalling onto a Schwarzschild black hole. The streamlines 
start at a spherical shell, where boundary conditions are fixed with wide flexibility, 
and are followed down to the point at which they either cross the black hole horizon 
or become incorporated into an equatorial thin disc. Analytic expressions for the 
streamlines and the velocity field are given, in terms of Jacobi elliptic functions, under 
the assumptions of stationarity and ballistic motion. A novel approach allows us to 
describe all of the possible types of orbit with a single formula. A simple numerical 
scheme is presented for calculating the density f ield. This mode l is th e relativistic 
generalisation of the Newtonian one developed bv iMendoza et ahl (l2009f ) and, due to 
its analytic nature, it can be useful in providing a benchmark for general relativistic 
hydrodynamical codes and for exploring the parameter space in applications involving 
accretion onto black holes when the approximations of steady state and ballistic motion 
are reasonable ones. 
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j^ ! 1 INTRODUCTION 

^— ^ . For several decades, accretion onto black holes has been re- 
cognised as the mechanism behind some of the most power- 
ful astrophysical phenomena (see e.g. iFrank et al.l l2002h . 
It has been extensively studied in the con text of Active 
Galactic Nuclei (AGN) dCenzel et al.ll2010l ). Gamma-Ra y 
Bursts (GRBs) (|Piranll2004l ), X-ray binari es |Kin3 [l99l) , 
compact binary coalescence (IHughesI 20031 and tid al dis- 
ruption of stars by black holes ( Rosswog et al.ll2009l ). 

Spherical accretion onto a black hole was found to have 
low efficiency for conv erting gravitat ional potential energy 
into emitted radiation (|Shapirdll973h and so rotation of the 
accreting matter has been invoked to give the increased ef- 
ficiencies required by observations. Rotation in the flow of- 
ten leads to a centrifugal hang-up with the formation of an 
accretion disc around the central object through which ma- 
terial inspirals slowly under the action of dissipative mech- 
anisms which can give efficient conversion of gravitational 
potential energy into radiation. The studies of accretion 
discs involve many physical inputs including magnetohydro- 
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dynamic turbulence, the behaviour of highly ionised and de- 
generate matter, radiative processes and radiation transfer, 
the formation of strong shocks, nuclear burning, etc; how- 
ever, gravity and rotation play the dominant role in determ- 
ining the overall accretion regime and bulk dynamics. 

Mo dern a c cretio n theory began with the pioneering 
work bv lBondil (|l952l '). in which he gave the analytical solu- 
tion for the steady spherical accretion of an ideal gas cloud 
onto a central object in Newtonian grav ity. The relati vistic 
extension of this was then developed bv lMichell (|l972h who 
took a Schwarzschild black hole as the central accretor. Ro- 
tating inflows were discussed by IPrendergast fc Burbidgd 
(196S '). who introduced the id ea of accretion discs which 
was then developed further b y IShakura fc SunvaevI (|l973) 
and lNovikov fc Thornel (|l973l ). 

Formation of an accretion disc by infall of a r otating 
gas cl oud onto a central object was flrst treated bv lUlrichI 
(| 19761 ) in the context of star formation and accretion discs 
around protostars. In that Newtonian analysis, the accret- 
ing gas was considered to start falling in from a spher- 
ical shell located infinitely far away from the central object 
where it had uniform density. Furthermore, it was assumed 
that the shell was rotating uniformly and that the fluid ele- 
ments were following parabolic trajectories. The disc forms 
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in the plane perpendicular to the cloud's rotation axis as 
a result of the collision of streamlines coming from mirror- 
symmetric points on opposite sides of the equator. (See e.g. 



Cassen fc Moosmanlll98ll:lLin fc Pringlelll990l : [stahler et alj 



19941 : iMendoza et alJIJOOi : lNagelll2007l . for a clear descrip- 
tion of Ulrich's model and the disc formation process.) The 
relativistic e xtension of Ulrich's model was in vestigated nu- 
merically by lBeloborodov fc lUarionovl (|200l|), in relation to 
wind-fed X- ray binaries, and h a s also been studied ana- 
Ivti callv by iHuerta fc Mendozal (|2007l ') (HM07 hereafter) 
and iMendoza fc Huertal ( 200S ). In the paper to which the 
present one is a follow-up ( Mendoza et al.ll2009l . - referred 
to as Paper I in the following), Ulrich's model was exten- 
ded, within Newtonian gravity, by considering a finite cloud 
radius and allowing for more general boundary conditions. 

In the present paper, we give a full relativistic general- 
isation of Paper I, studying the stationary inflow of a rotat- 
ing cloud of non-interacting particles around a Schwarzschild 
black hole. We shall assume that the particles follow ballistic 
trajectories and thus, that flow lines correspond to timelike 
geodesies of Schwarzschild spacetime. Building o n previous 
analytical studies (see e.g. IChandrasekhailll983l '). we intro- 
duce a novel approach which allows us to describe all of the 
different types of trajectory with a single analytic expres- 
sion. The model tracks the infall of particles injected from 
a finite spherical surface centred on the accretor where the 
density and velocity distribution of the particles are fixed 
with wide flexibility. Due to the rotation, a disc-like struc- 
ture forms in the equatorial plane. In the present work we 
focus on the accretion flow feeding the disc, and not on the 
behaviour of matter within the disc itself, where the ballistic 
approximation certainly breaks down. Our aim is to give an 
analytic description of this idealised accretion flow, starting 
far away from the central object and following it down to 
the vicinity of the equatorial disc, where we then stop our 
calculations. In the present work, the disc and the black hole 
are treated just as passive sinks of particles and energy. 

In this paper we follow an analytic approach similar to 
that in HM07, although in Section [7] we show that the fi- 
na l expressions of HM07 need to b e modified. Furthermore, 
in iBeloborodov fc lUarionovl (|200ll ') and HM07 the authors 
considered as the boundary condition a uniform-density, 
rigidly-rotating dust shell with all of the fluid elements fol- 
lowing parabolic-like motion. In contrast, the present model 
allows for more general distributions of density, rotation pro- 
file, accretion rate and particle energy. 

We propose to use this analytic calculation as a bench- 
mark to test the ability of general relativistic hydrodynam- 
ical codes in recovering geodesic motion when weakly inter- 
acting particles moving on a fixed background metric are 
considered. 

Moreover, the model presented here might also con- 
stitute a valuable tool for exploring the effect of different 
velocity and density distributions on the overall accretion 
scenario, in astrophysical applications where the ballistic 
and steady state conditions are reasonable approximations 
such as supersonic accretion from a w ind-fed X-ray binary, 
as in IBeloborodov fc lUarionovl (J200ll ); sub-Eddiri gton gas 



before full hydrodynamical simulations are performed, be- 
comes especially useful if one considers that in most of these 
situations the angular momentum distribution of the accret- 
ing matter is highly uncertain, even if it is clear that rotation 
does play a crucial role in determining the overall accretion 
efficiency. 

As an example of the validity of such an exercise, in 
Section [6] we compare the velocity and density fields, as 
predicted by the analytic model, with ones from a roughly 
equivalent numerical simulation from LR06. In that work, 
the authors explored numerically the accretion fiow of the 
collapsing interior of a massive star towards a newborn 
black hole, mimick i ng re lativistic effects by means of a 
IPaczviisky fc Wiital (|l980l ) pseudo-Newtonian potential. 

The structure of our paper is as follows. In Section [2] 
we present the model and its assumptions. In Section [S] a 
general expression for the fluid streamlines is given. In Sec- 
tion U the velocity fields are described as seen by observers 
located at infinity and by local observers. Using the continu- 
ity equation, a numerical scheme for calculating the density 
field is developed in Sectional The model is illustrated with 
a simple choice of boundary conditions and then compared 
against a numerical simulation from LR06 in Section |6| In 
Section[71 we adopt instead the boundary conditions of the 
lUlrichI (|l976l ) model, and give its general relativistic exten- 
sion. Finally, in Section [S] the non- relativistic limit is con- 
sidered and the results given in Paper I are recovered. 



2 MODEL DESCRIPTION 

We are interested here in modelling a rotating cloud of 
particles falling towards a central black hole with mass 
Af, whose exterior gravitational field is described by the 
Schwarzschild metric. We assume that the accretion flow has 
reached a stationary situation characterised by a constant 
accretion rate M, where from now on, a dot denotes a de- 
rivative with respect to the proper time r. Additionally, we 
assume that the gravitational fleld of the black hole is the 
main factor determining the fluid dynamics. We therefore 
neglect the effects of ffuid self gravity, pressure gradients, 
fluid viscosity, radiation pressure, neutrino capture, etc. In 
other words, we give a ballistic treatment of the fluid flow. 
Consider a cloud subdivided into equal mass fluid ele- 
ments. We let Q be the matter density and n be the ba- 
ryon number density, and introduce an average baryonic rest 
mass, (m), such that g — (m) n. We assume that fluid ele- 
ments are continuously injected from a spherical shell of 
radius Tq. This shell represents the outer boundary of the 
model where the fluid properties are fixed as 



Qo^ Q {ro, do, <j)o) , 
ro = r (ro, 9o, 4>o) , 
Oo = (r-Q, ^0, ^o) , 
00 = [ro, do, (t>o) , 



(2.1) 
(2.2) 
(2.3) 
(2.4) 



accretion onto galactic nuclei (see e.g. Blae g 12007*): or col- 
lapsing stellar cores , as in Lee fc Ramirez-Ru iz (200 6|) (LR06 
hereafter) (see also lLopez-Camara et al.ll2009l : ITavlor et al.l 
[201l|). Exploring the role of different boundary conditions. 



with r, 6 and (j> being the radial, polar and azimuthal ve- 
locity components, respectively. Figure [l] shows a schematic 
diagram of the accretion scenario. 

We take the four distribution functions in eqs. (|2.ip - 
(|2.4|l to be differentiable and to be symmetric with respect 
the equatorial plane, i.e. 
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Figure 1. Schematic illustration of tiie present model. Boundary 
conditions are fixed at a shell with radius ro and fluid elements 
then move along ballistic trajectories. We assume an axisymmet- 
ric distribution of angular momentum with its magnitude increas- 
ing monotonically towards the equator. The streamlines can be 
divided into three groups: (i) ones with low angular momentum 
that accrete directly into the black hole; (ii) ones with interme- 
diate angular momentum which arrive at the equator but do not 
find a circularisation radius; (iii) ones with large enough angular 
momentum to be incorporated into an equatorial Keplerian-type 
accretion disc. 



/(eo) = /(^-eo), 



(2.5) 



with / = Qo, r-o, 00, (pg. 

Using the boundary distributions for ro and qq, one cal- 
culates the total accretion rate as 



M = 



-ro 



ZTT TT 

// 



QoTo sin6'od6'od(?!>o. 



(2.6) 



Given the ballistic approximation, we have that the spe- 
cific energy, E, and the specific angular momentum, h, are 
conserved along streamlines. Clearly their actual values vary 
from line to line and will be distributed in a certain fash- 
ion (see eqs. 13.211 and I3.22|) as a function of the boundary 
conditions in eqs. (|2.1|l - (|2.4p . 

It is a feature of the Schwarzschild spacetime that, for a 
given mass of the central object and a given streamline, with 
energy E{0o, 4>o), there exists a critical value of the specific 
angular momentum|j /ic, such that if /i(^o, 4>o) < he then the 
corresponding fluid line crosses the black hole horizon (loc- 
ated at Ts = 2 GM/(?) before reaching the equatorial plane. 
On the other hand, for h{6o, 4>o) > he the streamline reaches 
the equator and encounters there an equivalent streamline 
coming from (tt — 6o, <j)o)- In a real physical situation one 
expects that if these collide supersonically then two shock 
fronts will appear above and below the equator, with the 



Note that taking he = 2rsC as the critical value, is valid only 
for streamlines with parabolic-like energies, i.e. E = (? . 



fluid being incorporated into a disc-like structure. Provided 
that there is an efficient dissipation mechanism, the shock 
fronts will remain pinned down to the equator and the disc 
will remain thin. It is clear that the study of the disc dy- 
namics requires a full hydrodynamical treatment, in which 
redistribution of angular momentum and energy losses are 
self consistently taken into account, but such an analysis 
lies beyond the scope of the present work. Instead, we shall 
just assume here that an efficient mechanism dissipates all 
of the kinetic energy associated with the vertical component 
of the velocity, in such a way that an infinitesimally thin disc 
forms in the equatorial pl ane which is then taken to act as 
a passive energy sink. See iLopez-Camara et al.l (|2009r ) for a 
full hydrodynamical simulation of a coUapsar in which they 
show that an isothermal disc does indeed remain thin. 

In principle one could relax the condition in eq. (|2.5|l 
and not assume any particular symmetry for the fluid prop- 
erties at the boundary. However, in that case we would not 
have the symmetric collision of streamlines described above. 
This would lead to formation of a warped disc, making the 
situation much more complicated. 

Following LR06, the streamlines can be divided into 
three groups depending on the value of their specific angu- 
lar momentum (see Figure [ij : (i) streamlines with low an- 
gular momentum, which go directly into the black hole; (ii) 
streamlines with intermediate angular momentum, which 
form a small disc within which accretion proceeds on a free- 
fall time scale ( this situation corresponds to the one explored 
numerically bv lBeloborodov fc Illarionovll200ll ') ; (iii) stream- 
lines with larger angular momentum, which have sufficient 
centrifugal support on their arrival at the equator so that 
subsequent accretion would occur on a viscous time scale in 
a Keplerian-type accretion disc. 

Now consider the situation in which two neighbouring 
streamlines start approaching each other in such a way that 
they would intersect. This type of encounter is qualitatively 
different from the head-on collision described above, since 
here, and with a full hydrodynamical treatment, the two ap- 
proaching streamlines would be prevented from intersecting 
by the smooth action of pressure forces. It is clear, however, 
that this cannot be handled within the ballistic approxim- 
ation and so we must restrict our analysis to distribution 
functions for which streamlines do not cross before reaching 
the equator. 

Provided that there are no early intersections, then we 
can use the initial angular position, (Q^, <po), as a label of 
the individual streamlines. In the next Section we give an 
analytic expression for the streamlines that, for any given 
radius r, constitutes a mapping {Og, (f>o) -f-^- (9, 0). The con- 
dition of no intersection is equivalent to requiring that the 
Jacobian of this transformation should be non-singular, i.e. 



'-^M 



d(t> 



de 



d<t) 



>0. 



(2.7) 



While formal, the actual evaluation of this for given bound- 
ary conditions is not a trivial task in general. However, if 
there is axisymmetry it is a sufficient condition that the 
specific angular momentum should be a non-decreasing func- 
tion of 9o going from the polar axis to the equator, and that 
its magnitude should be sufficiently small so that no fluid 
elements reach a turning point in their trajectories before 
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crossing the equator. For more general cases, it seems best 
to proceed by trial and error. 



3 STREAMLINES 

Under the ballistic approximation, the streamlines of the ac- 
cretion model correspond to trajectories of freely-falling test 
particles (i.e. timelike geodesies) in a Schwarzschild space- 
time. In this Section we give a general analytic expression 
for the spatial projection of these geodesic lines. 



3.1 Description in the orbital plane 

Let us consider a given fluid element starting to fall in from 
the point {ro, 6a, (j>o) ■ We denote by O the general frame of 
reference with coordinates (f , r, 9, (j>) having its polar axis 
coinciding with the rotation axis of the cloud. 

In a general spacetime, the geodesic equations consist of 
a second order, non-linear, coupled system of four differen- 
tial equations. Nevertheless, with the aid of the underlying 
symmetries of the Schwarzschild spacetime, it is found that 
the corresponding geodesic equations become separable and 
red uce to a set of four first order differential equations (see 
e.g. iFrolov fc Novikovlll998l ). A further simplification com- 
ing from the spherical symmetry is that the particle motion 
is confined to a single plane. We denote by O' a frame of ref- 
erence with coordinates {t,r,-d,ip), especially adapted such 
that its equatorial plane (i? — 7r/2) coincides with the orbital 
plane. 

From the view point of O' , the geodesic equation asso- 
ciated with the polar angle becomes i? = 0. This confirms 
that the whole trajectory stays in a single plane. On the 
other hand, the equations corresponding to the time and 
azimuthal coordinates lead to two first integrals of motion, 
namely the relativistic total specific energy, E, and the spe- 
cific angular momentum, h, defined as; 



E- 



2 . 

: r if, 

r 



(3.1) 
(3.2) 



With the aid of eqs. (|3.1|) and (|3.2p together with the 
normalisation condition for the 4-velocitjo (u'' = dx^/dr), 
u^u^ = —c?, one gets the following equation governing the 
proper time evolution of r for the fiuid element 



= e + 



2GM 






where we have introduced the re-scaled energy 



(3.3) 



(3.4) 



This definition of energy is convenient because it facilit- 
ates comparison with the Newtonian case since, in the non- 
relativistic limit, e converges to twice the Newtonian total 
energy. 

Since we are interested in a stationary regime, it is con- 
venient to recast the time derivative in eq. (|3.3p as a derivat- 
ive with respect to the azimuthal angle ip (d/dr = ipA/Aip). 
Doing this together with the aid of eq. 1)3. l[l allows us to 
rewrite eq. (|3.3p as 



Ar 

dip 



VW) 



with 



7^(r) 



2 GMr 



r)^ 



(3.5) 



(3.6) 



The minus sign in eq. (|3.5|l is needed because, in the current 
accretion scenario, ifi increases while r decreases as the fiuid 
elements approach the equator. 

Being a fourth degree polynomial, TZ{r) has four roots 
(with possible multiplicity) and we can write it as 



TZ(r) = e(r — ra){r — rb){r — r^){r — r^)- 



(3.7) 



Explicit expressions for the roots are given in the Appendix. 
It is clear that r = is one of them; we treat it in the 
same way as the others though, since keeping an explicit 
reference to it allows us to give a general expression for the 
streamlines. Given that one root is zero, then necessarily at 
least one other root must be real. The remaining two roots 
can either also both be real, or can form a complex conjugate 
pair. 

After some analy sis of T^{r) (see, e.g. 
iMielnik fc Plebanskil Il962l i it follows that if £ < 0, all 
of its real roots are non-negative, while for e > 0, there 
is exactly one negative root. For e = 0, one of the roots 
has diverged to infinity and TZ{r) reduces to a third 
order polynomial. This case is treated in further detail in 
Section [T] 

From eq. (|3.5|) it is clear that the radial motion is restric- 
ted by the condition TZ{r) > and so r is either bracketed 
in between two consecutive real root^ (bounded motion) or 
is unbounded above. 

Let Va be the lower bound (periastron). (This includes 
the possibility r„ — 0, which corresponds to a plunge orbit.) 
If the four roots are real and the motion is upper 
bounded, we let r^ be the upper bound (apastron), while 
for unbounded motion, we let r,, be the only negative root. 
We then take the two remaining roots, r,, and r^, such that 
r-c < Ta- 
li two of the roots are complex, we call those roots 
Ti, and r^, with r^ = r^*, and let r^ be the remaining real 
root. This case, with complex roots, represents the purely- 
relativistic "capture" type of orbit, for which particles 
with non-vanishing angular momentum do not have strong 



^ Throughout this paper we use Einstein's summation conven- 
tion with Greek indices running over spacetime coordinates while 
Latin indices run just over the spatial components. 



•^ Whenever these bracketing roots are positive, they correspond 
to turning points for which the radial velocity vanishes and the 
particle goes from moving inwards to moving outwards, or vice 
versa. 
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enough centrifugal support to prevent them from faUing dir- 
ectly into the black hole. 

We set the origin of ip at the periastron of the orbit, i.e., 
r{(p = 0) = r„. The solution of eq. (|3.5p is then equivalent to 
finding solutions to the following quadrature problem: 



dr' 



VW 



h' 



(3.8) 



This integral is solvable in terms of Jacobi elliptic functions, 
as has been discussed several t i mes in the literature (see , 



g.lHagiharalll93(]l:lDarwinll 19591: iMielnik fc Plebanskill 19621 : 



lMetznerlll963l : IChandrasekhaiJll983l : iMiro Rodriguez! ligSTl ). 
These works have followed several different approaches, fo- 
cusing on different aspects of the solution and using different 
notations. Labelling the roots in the way described above en- 
ables us to give a novel description of the different types of 
trajectory by means of a single analytical expression, which 
summarises the results found in the previous work: 



n(rd - Ta) - ra{rb - r„)cn^(^, fc) 



rd- 



{rb — ra)cn2(^, fc) 



with 



^ 



•^/e(r„- r^){rd- n) 
2h 



■V', 



(3.9) 



(3.10) 



where en (^, fc) is a Jacobi elliptic function with modulus fc 
given by 



fc = 



(u- r^){ra- Tc) 
{Td- rb){r^ -Ta)' 



(3.11) 



Note that, though general, the expressions in eqs. (|3.9|l - 
(|3.11|l should be handled with care in two particular cases: 
when complex roots are involved and when e — >■ 0. In the first 
case, some intermediate terms will be complex quantities 
even though the final result will always be a real number. In 
the Appendix we give alternative expressions for this first 
case, while the second case is discussed in Section [T] 



the one in O by an angle 7r/2 — 9a and performing a series of 
rotation operations|j one obtains the following relationships 



COs(</3 — ipa) 



COS 6 

COS 9a 

cot 9 
cot 9a ' 



(3.13) 
(3.14) 



Differentiation of eqs. (|3.13p and 1)3. 14[) with respect to 
the proper time gives 



sine* 



x/cos^ 9a — cos^ 9 

. ■ sm9a9 

sin9cl>=—= (3.16) 

Vcos^ 9a — cos^ 9 

Evaluating eq. (|3.16[l at ^o gives the following expression for 
9a in terms of known quantities 



sm tJa<po 



?o2-Hsin^0oV 



(3.17) 



From this last equation it is clear that if 9q = then 9a ~ 9o. 
On the other hand, evaluation of eqs. (|3.13|) and (|3.14|) 
at {9o, 4>o) allows us to calculate ipa and cj>a through 



Va = fa — COS 



= ©0 — COS 



COS) 



COS^a 

1 / COtgp 

cot6l„ 



(3.18) 
(3.19) 



where, from eqs. 1)3. 91) and p.lO|) . 
-2/i 



fo - 



\/e(ra -~ra){rd- u) 



{r-d- ra){rb-ro) 
{rb~ ra){rd-'ro)' 



(3.20) 



Using eq. (|3.1|) and combining eqs. p.l5|) and p.l6p . we can 
express h in terms of the boundary conditions as 



3.2 Relation between the frames of reference 

In this subsection we relate the descriptions made in O and 
in O' , no ting that, for fixed t and r, the Schwarzschild metric 
(see, e.g. iMisner et al.lll973l l reduces to that of an ordinary 
3-sphere and, hence, basic rotation operations can be per- 
formed. The easiest way to relate ip to the angles 9 and 
is by means of introducing the turning point in the polar 
motion, i.e., a point 9a for which 



.) = 0. 



(3.12) 



The polar motion of a non-equatorial trajectory is always 
characterised by two turning points which are symmetric 
with respect to the equator. That is, if 9a satisfies eq. (|3.12p 
then IT ~ 9a does as well. Here we choose 9a in the same 
hemisphere as 9o- 

Noting that the polar axis in O' is tilted with respect to 



h^ri\l9,? + smH„4>o'', (3.21) 

while evaluating eq. (|3.3[) at r = Tq allows us to calculate e: 



.2 2GM h 

e = ra 1- ^ 

To ro-= To 



rsh_ 

3 



— -^. (3.22) 



From eqs. (|3.2ip and (|3.22|l it is clear that h and e are 
functions only of the initial angular position {9o, (po) (i.e. 
they depend only on the boundary conditions) and that they 
are conserved along fluid lines. 

The expressions given in this Section constitute a de- 
scription of the fluid lines as a function of the relat- 
ively simple but still general boundary conditions given in 

■* First make a rotation of —<fia about the z axis in the C frame, 
then one of tt/2 — 9a about the resulting y axis, followed by one 
of (l>a about the final z axis. 
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eqs. (|2.ip - (|2.4|) . What one needs to do to use these expres- 
sions in practice is as follows; from the boundary conditions 
find the values of 6a,, h and e using eqs. (|3.17p . p.2ip and 
(|3.22|l . respectively. This determines the coefficients of the 
polynomial TZ{r), as given by eg. (|3.6|l . and its roots can 
then be calculated in the way described in the Appendix. 
The next step is to find the numerical values of k, ipo, ip^ 
and (j>a and finally one can make use of eqs. (|3.9|) . (|3.10|) and 
(|3.13|) to track a given fluid element in its descent from Tq 
towards the equatorial plane, i.e., as 9 sweeps from 9o to 
7r/2. 



4 VELOCITY FIELD 



Using eqs. ((3TT1) - (p^ together with eqs. (pT5|) and (|3TT6)) . 
one obtains the following expressions for the velocity field 



u = — ■ 1 



e + 



2 GM K^ 



('-t). 



u = ± 



h\/cos^ da — cos^ 9 



h sin6'„ 



(4.1) 

(4.2) 

(4.3) 
(4.4) 



where the sign in eq. (4.3) is positive for < ^ < 7r/2 and 
negative for 7r/2 < 9 < n. 

It is convenient to introduce as well the velocity field 
as described by a set of local observers, since their measure- 
ments correspond to physical intervals of proper distance 
and proper time. Associated with each local observer one 
has an orthogonal tetrad of 4-vectors constituting a local 
Minkowskian coordinate set of basis vectors. Let us denote 
with a bar, coordinates in the local frame of reference, i.e. 
let u^ be the components of the 4- velocity of a fluid element 
passing by a local observer; 



i^ = [jc,'yVj, 



(4.5) 



where 7 is the general relativistic Lorentz factor between the 
local observer and the fluid element, defined as 



1- 



V-V 



(4.6) 



E_ 

72 



('-?) 



and where 



r 


c\^ 




E ' 


e 


r u 




7 


4> 


r sin 9 u'l' 



7 



are the components of the three- velocity V . 



(4.7) 
(4.8) 
(4.9) 




Figure 2. Schematic illustration of the construction of a regular 
grid in {Oq, r) for calculating the density field. For clarity, the 
azimuthal direction is not included. The highlighted region con- 
stitutes one of the streamline tubes involved in the derivation of 
eq. 15211 • 



5 DENSITY FIELD 

The expressions for the streamlines and the velocity field 
given in Sections |3] and |4] are independent of the value of 
the density at the boundary and hence the scale for the 
density can be set arbitrarily. In this Section we derive a 
numerical scheme for calculating the density field based on 
the continuity equation. For a general curved spacetime, this 
equation can be written as 



V^(nM") = 



d 



—g dx^' 



^nu") =0, 



(5.1) 



where g — det[(jf,„] is the metric determinant, given in the 



Schwarzschild case by ; 



-r sin 9. Using the stationarity 



condition, eq. (|5.1|) reduces to 



d 



r^ sin 9 dx 



-(r sm9nu') = 0, 



(5.2) 



which is simply saying that the spatial divergence of the 
particle number fiux 3-vector nu^ is zero. We can integrate 
this equation over a streamline tube, i.e., a volume element 
consisting of a collection of streamlines coming from an area 
element da|ro at the shell rg and ending up at a second sphere 
with arbitrary radius r < Tq. Such a streamline tube is illus- 
trated in Figure [2I By means of Gauss's theorem it follows 
that 



nu'daA = nu'dai 



(5.3) 



The differential area element orthogonal to the radial 
direction is given by 

da,. = r^sin6ld0d(;/), (5.4) 

and so eq. 1)5. 3|l can be written as 

no UoVg sin 9o d9o d4}o — nu^ r sinSd^dtji. (5.5) 

Using the relation d9d(j) = Jd9od(j}o, where the Jacobian J 
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was given in eq. (|2.7p , we get the required result for the dens- 
ity field 

no -Uo To sin So ,_ „, 

As long as J > 0, which means that no intersections 
of the streamlines occur, and w' < 0, which implies that 
no turning points in the radial motion exist up until the 
equator is reached, the expression for the density field given 
by eq. (|5.6p is well defined and has no singularities. 

The partial derivatives required in the calculation of 
J represent a complex computation involving derivatives of 
an elliptic integral with respect its argument, modulus and 
integration limit. On the other hand, it is straightforward 
to evaluate it numerically and so it does not seem worth 
searching further for a full analytic expression. 

We construct a suitable grid for calculating J in the fol- 
lowing way. We start with a homogeneous partition of the 
initial angles (^o, 4>o) and then follow the fluid lines down 
to the equator in regular radial steps. At every grid point 
(9o, 00, t) we store the values of 9, (p and u^ and then cal- 
culate J by means of standard finite differences. Figure [2] 
illustrates the construction of such a grid. 



6 EXAMPLE MODEL 

Here we illustrate the accretion model presented in the pre- 
vious sections by applying it with the boundary conditions 
considered in Paper I, i.e. ones for a uniform shell of matter 
in uniform rotation: 



Qo = const, 
f-Q = const, 
(j)o ~ const, 

^0 = 0. 



(6.1) 
(6.2) 
(6.3) 
(6.4) 



This choice leads to several simplifications, the most import- 
ant being that the accretion flow is then axisymmetric and 
so, the Jacobian of the angular transformation in eq. (|2.7|) 
simplifies to 



J 



dOo' 



(6.5) 



From eq. (|2.6p . we have that the total accretion rate is given 

by 



M = 4Tvr„^eo\ro\, (6.6) 

and, by substituting the boundary conditions in eqs. (|6.ip - 
(|6.4p into eq. (|3.2ip , we get the following distribution of spe- 
cific angular momentum: 



^(^o) = fee sin ^0, 



(6.7) 



where h^ — Tq^q is the maximum value of the specific an- 
gular momentum, which is reached at the equator of the 
shell. Since 6a = 0, we have that for every streamline 
da = do, 'Pa = Pa and (pa = (po givlug a simplification of 
eqs. (|3.13p - p.l6p . The velocity field is again described by 



the expressions given in Section [4] after making the substi- 
tution 9a = ^0- Regarding the density field, from eqs. (|5.6p 
and (|6.5p . we find the expression 



no Mq ^0 sin 6a { 86 



W r^ smt 



(6.8) 



Figure [3] shows the projection onto the R-z plane {R = 
r sin 6, z = r cos 6) of the streamlines, velocity field and dens- 
ity contours for three different combinations of f-g and h^. In 
each case we have set ra = 20 Ts . 

For the choice of boundary conditions being used here, 
the radius of the outer edge of the disc formed as matter 
reaches the equatorial plane, Vn, can be calculated from 
eq. ((3T0)) and (f3T3l) . taking first 6 = tt/2 and then 6o = 7r/2, 
giving 



c« = 



^£(r-„-re)(rd-r6) 
2K 



f ^ 

[po+- 



(6.9) 



and then substituting the result into eq. (|3.9p . In Figure |3J 
we have plotted ro as a function of /i, as obtained from 
eq. (|6.9p for a fixed value of ro and five different values of 
To- Note that ro is the radius at which the particle with the 
largest angular momentum first impacts on the equatorial 
plane; no further motion is then considered. For sufficiently 
small rjj, the "disc" will be one within which accretion into 
the black hole then proceeds on a free-fall time scale, as 
discussed earlier, and only for larger values of ro could it 
give rise to a Keplerian-type disc. Naturally, ro increases 
monotonically with increasing ft,, while it decreases with 
increasing Iroj. For a given specific angular momentum, ro 
can be substantially smaller when larger values are taken for 
I'^ol. 

Let us now consider boundary conditions corresponding 
to the collapse of a massive stellar core such as those stud- 
ied numerically by LR06. In that work, the authors investig- 
ated the formation of an inviscid, small-scale accreting disc 
by making a 2D Smoothed Particle Hydrodynamics (SPH) 
simulation with rel a tivist ic effects being mimicked with a 
IPaczviiskv fc Wiita (|l980t ) (PW) pseudo-Newtonian poten- 
tial. They included a detailed treatment of neutrino cool- 
ing and considered an equation of state with contributions 
from radiation, e pairs, a particles and free nucleons. We 
have made comparisons for several of the models discussed 
in LR06, finding similar results in all cases. For illustration, 
in Figure[5]we show just one of them, namely the numerical 
run in which LR06 took a central black hole of 4Mq and an 
external spherical boundary at ro — 50 r^ from which SPH 
particles were continuously injected with a constant accre- 
tion rate of M = 0.01A/q/s. As radial infall velocity, they 
took the velocity of free fall from infinity, i.e. To = — \/l/50 c, 
and the specific angular momentum of the fluid elements at 
ro was assumed to follow a rigid body rotation distribution 
with a maximum of h^ = 1.9 rsC at the equator of the shell. 

The top left panel in Figure [5] shows the accretion flow 
as calculated from the analytic model, while the top right 
panel shows a late-time snapshot of the LR06 simulation, 
taken when a quasi-stationary state had been reached. In 
general, there is quite good agreement between them until 
the streamlines have approached the equatorial disc. This 
is not surprising, bearing in mind that the flow is highly 
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i-o = -0.001 c , /i, = 2 CjC 





g 


m 


m 






~~''~ \ j^^ 0.648/c- 



n 


wii 



-20 



-10 







10 



20 



Figure 3. Streamlines, velocity field and density contours for 
different flow parameters. In each of the panels we have taken ro = 
20 rs , while the values of the remaining parameters are indicated 
in each plot. The plots are projections onto the R-z plane and 
the colour coding corresponds to the logarithm of the particle 
number density, Log(n/no). The arrows correspond to the V" 
and V components of the velocity field. 




/ijr,c 



Figure 4. The radius of the outer edge of the disc formed as 
matter reaches the equatorial plane, rn, is plotted against the 
specific angular momentum at the equator, h^. A fixed value 
of ro = 20 Ts is considered. From top to bottom, each curve 
corresponds to a different inward radial velocity with values 
\ro\/c = 0, 0.2, 0.4, 0.6, 0.8 respectively. 



supersonic all the way down to the vicinity of the equat- 
orial disc, where strong shocks then appear in the hydro- 
dynamical simulation because of collisions between stream- 
lines coming from opposite points in the two hemispheres. 
The other four panels present a detailed comparison of the 
spatial components of the velocity and the density at four 
spherical cuts. Here, we see very good agreement between 
the analytical and numerical results for u® and u*^. For it'' 
and g, there is quite good qualitative agreement, although 
the numerical results for q suffer from numerical noise in- 
herent in the interpolation scheme at low particle number 
densities, and the numerical results for u^ show higher radial 
infall velocities, as expected given the use of the PW pseudo- 
Newtonian potential there, which artificially enhances the 
radial acceleration. 



7 RELATIVISTIC EXTENSION OF ULRICH'S 
MODEL 

In this Section we adopt bound a ry conditions correspond- 
ing to the models of lUlrichI (|l976l ). lBeloborodov fc lUarionovl 
(|200H ) and HM07, i.e. we consider infall from an infinitely 
large spherical shell of matter with all of the fiuid elements 
having parabolic-like energies (ro — >■ oo, ro = and hence 
e = 0). 

In Ulrich's model, the location of the outer edge of the 
equatorial disc is simply related to h^ through 



rk = 



GM' 



(7.1) 



In Newtonian gravity, r^ corresponds to the Keplerian radius 
of a circular orbit with specific angular momentum h^ as 
well as to the semi-latus rectum of a parabolic orbit with 
this specific angular momentum. 

Note that a vanishing initial radial velocity implies r^ = 
To and, consequently, re — while 



h' 



ra.d 



AGM 



1±a/1-4 



\ir) 



(7.2) 
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Figure 5. Comparison between the analytic model and one of the LR06 simulations. The plots are for an accretion flow towards a 
black hole with mass M = 4Mq , starting from a spherical shell at ro = 50 r^ where the matter density and radial inward velocity are 
uniform and given by go = 5.29 X 10^ g/cm'' and r-Q = —-^1/50 = —0.14 c, respectively. The specific angular momentum distribution 
at the shell corresponds to uniform rotation with a maximum of h^ = 1.9 rsC. With these boundary conditions, the total accretion rate 
is iW = 0.01Mq/s. The top panels show a projection of the accretion flow onto the R-z plane, together with isodensity contours of 
the analytic solution (left) and the LR06 numerical simulation (right). The remaining four panels show the velocity components and 
the density at the radial cuts r/vs = 20, 15, 10, 5 with the analytic and numerical results being represented by solid and dashed lines 
respectively. 
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Figure 6. Streamlines for tlie rclativistic extension of Ulrieh's 
model [ro — > oo, ro = 0). The continuous line corresponds to 
ht; '^ TsC, i.e. to the non-relativistic limit in which Ulrieh's model 
is recovered. The dotted line corresponds to h^ = TsC while the 
dashed line is for h^ = he ~ 0.754 TsC, in which case the disc forms 
entirely inside r^, i.e. ro = Ts- The quarter-circles correspond the 
location of the event horizon in each case. 



Taking the limit ro — >■ oo in eq. (|3.9|) - (|3.1ip and noticing 
that £ ro — >■ —2 GM, one gets the following expression for the 
streamlines 



-rdsn^(^, k) 
;n2(^, k) ' 



with 



c = 



GMr^ ifi 



(7.3) 



(7.4) 
(7.5) 



Figure [6] shows plots of the streamlines for different val- 
ues of fee. From these, we see how the radius of the equat- 
orial disc decreases as h^ decreases starting from ro — r^ 
for hc^ rsC (corresponding to the disc radius in the Ulrich 
model) down to ro = rs when h^ — hc~ 0.754 rsC. 

HM07 followed an analytic approach similar to the one 
presented here, but used an incorrect transformation for the 
polar and azimuthal angles between the systems of reference 
O and O' . In particular, eq. (27) in HM07 is not correct and 
should be substituted by eqs. (|3.13|l and (|3.14|) of the present 
paper. 



8 NON-RELATIVISTIC LIMIT 

In this Section we consider the non-relativistic limit, for 
which /le 3> TsC. By taking this limit in eq. (|3.22p . it follows 
that e becomes equal to twice the total Newtonian specific 
energy: 



■24. '^ 



2GM 

7 

ro 



(8.1) 



and, from eq 
7^(r) (i.e., r, = r^ 
real and given by 



one has that r = is a double root of 
- 0), with the remaining two roots being 



GM 



(-lie). 



(8.2) 



where e = ^l + £{h/GMy. 

From eq. p.llf) it follows that fc = 0. For a null value of 
the modulus, the Jacobi elliptic functions reduce to ordinary 
trigonometric functions: 



cn(a;, 0) = cos(a:), sn(a;, 0) = sin(a::) 

and so eqs. (|3.9|) and (|3.10|) become 



■3) 



TaT-b 



"a + {rt — ra) cos'^{tp/2) 
h^/GM 



1 + e cos if 



(8.4) 



The second equality in eq. (|8.4|l is the well-known expres- 
sion for motion along a conic section with eccentricity e, 
representing all of the possible types of orbit in Newtonian 
gravity. 

All of the expressions derived in Section [3.21 concerning 
the angles tp, ipo, 'fia, 0, Oo and 9^ remain valid, but note that 
the expression for the initial orbital angle is simplified: 



ifo ~ COS 



h^/GM - ro 



ero 



(8.5) 



The velocity field is given by eqs. (|4.2p - (|4.4|l after taking 
Vs = and noticing that, within the non-relativistic limit, 
proper time intervals become identical with coordinate time 
intervals: dr = di. The Newtonian results presented in Pa- 
per I can be recovered from the expressions given here once 
the boundary conditions in Section [6] have been adoptedlj 
The density field is again given by eq. (|6.8p . Note that for the 
particular boundary conditions considered in Paper I, dir- 
ect evaluation of J is straightforward and gives the analytic 
expression presented there. 



9 SUMMARY AND DISCUSSION 

In this paper we have presented an analytic solution for the 
streamlines of pressureless matter being steadily accreted 
towards a Schwarzschild black hole. The accretion flow is 
taken to start from a spherical surface far away from the 
central mass, with a wide range of boundary conditions be- 
ing used. The fluid streamlines are then tracked down to 
the point at which they either become incorporated into a 
thin equatorial disc or pass inside the black hole event ho- 
rizon. We have presented a simple numerical algorithm for 
calculating the density field. 

In our model, the fiuid streamlines correspond to time- 
like geodesies of Schwarzschild spacetime. Using Jacobi el- 
liptic functions and some standard identities, we have de- 
veloped a novel approach for describing all of the different 
types of trajectory with a single analytical expression. 



^ Note that in Paper I, the azimuthal angle was measured starting 
from the apastron instead of from the pcriastron as we do here. 
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This work extends the Newtonian model presen- 
ted in Paper I to a general relativistic one in Schwar- 
zschild spacetime and constitutes the analytic solution 
corresponding to the s c enario studied numerically by 
iBeloborodov fc Illarionovl (|200ll '). We are currently working 
on the further generalisation of the model to Kerr spacetime. 
For appropriate bound ary conditions , the present solution 
generalises the classical lUlrichI (|l976l ) model and gives cor- 
rected expressions for the HM07 model. We have shown that 
our model recovers the well-known Newtonian expressions in 
the non-relativistic limit. 

Our analytic solution can be used as a benchmark for 
testing general relativistic hydrodynamical codes. Clearly, 
such codes should be able to correctly reproduce geodesic 
motion for a cloud of non-interacting particles in a fixed 
metric. 

Despite the fact that the present model leaves out many 
physical processes which are relevant for the study of real- 
istic accretion onto black holes, it does include the two main 
factors determining the bulk dynamics: gravity and rotation. 
Furthermore, given its flexibility for setting boundary condi- 
tions, it can be useful as a computationally inexpensive and 
efficient tool for exploring the parameter space in applica- 
tions where the ballistic and steady state hypothesis are ap- 
proximately valid, such as: sub-Eddington accretion towards 
a supermassive black hole in a galactic nucleus, wind-fed 
and Roche-lobe fed X-ray binaries, and collapsar GRB pro- 
genitors. In this way one can gain physical insight and get 
order-of-magnitude estimates of the energy budget before 
undertaking full-scale simulations. Having such an explorat- 
ory tool can be especially relevant if one bears in mind that 
the parameter domain for many of these systems is vast and 
often uncertain. An example is the comparison made in Sec- 
tion E] with the LR06 SPH simulation, where rather good 
agreement was found between the analytic model and the 
numerical results. 
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APPENDIX : SOLUTION OF THE RADIAL 

EQUATION 

In this appendix we discuss the solution for tlie radial mo- 
tion in more detail. Explicit expressions for the roots of the 
fourth degree polynomial TZ{r) defined in eq. (|3.6|) are easily 
given in terms of the followi ng auxiliary quantities (see e.g. 
lAbramowitz fc Stegunlll970l ) 



= (2GAf)^ + 3e/l^ 



R = {2GMY + 9er ( GM + ^r, e 



D' 



R^ - Q\ 



(A.l) 
(A.2) 
(A.3) 



When D < 0, all of the roots are real and are given by 



Ti ■ 



r2 = 



rs = 



) cos I * ) - GM 



57 [^ 

s[^'="(*+f)-H' 

9 r 

''4= ^ \VQcOs{^ +TT 

with '5 being defined through 



GM 



cos(3*) 



R 

13/2- 



(A.4) 
(A.5) 

(A.6) 

(A.7) 

(A.8) 



When r4 < 0, we interchange r2 and r^ in order to satisfy 
^2 < Vs. In this way, when real, the four roots are ordered as 



In this way the solution to eq. p.Sp is always given by 
eq. (O. 

On the other hand, in the case D > 0, we take 



ra^ri, rb = r2, r^ = rg, r^ = r4, (A. 17) 

and once again the expression in eq. (|3.9p is a formal solu- 
tion to eq. (|3.8p . However, direct evaluation of this involves 
the use of complex quantities as intermediate steps. It is 
possible to rewrite eq. (|3.9[l as an expression involving just 
real quantities. For doing this, we introduce the following 
two real constants (having in mind the fact that ri = 0) 



a = ±^y{r^-r2)ir^~rs), (A.18) 

P = ^Av^, (A.19) 

where the sign of a coincides with that of e. We then define 



C2 



h 



^ = 2V^^, 



,2_ {a + pr-r^ _ {iTkr 



4q/3 



=F4fc 



(A.20) 
(A.21) 



where the sign accompanying k in last equations is the op- 
posite of the one for a . If we now invo k the following identity 
for Jacobi elliptic functions (see, e.g. lHancocklll917^ 



2 1 

=p fcsn (^, fc) 



en (^2, ki) 
1 + en (^2, fc2) ' 

we can rewrite equation (|3.9[l as 



(A.22) 



Ti < r2 < r3 < r4 



(A.9) 



On the other hand, for D^ > two of the roots form a 
complex conjugate pair. Making the definitions 



S^^/D-R, 
T = ^D + R, 

the non-zero roots in this case are given by 



(A.IO) 
(All) 



^^ /3r4[l-cn(C2, fca)] 

/3 - a - (a -f /3)cn (^2, fc2) ' 

Note that, in this case. 



Vo = 



-h _i 
: cn 



feap 



Pri + {a- P)ro 
Pn-{a + P)ro' 



(A.23) 



(A.24) 



In summary, the expression for the trajectory of a time- 
like geodesic in Schwarzschild spacetime can be written as 
follows: 



r2, ,3 = -^ [t - S - 4GM ± iV3{S + T) 



U=—{S-T-2GM), 



(A.12) 
(A.13) 



Since the radial motion is constrained to satisfy TZ{r) > 
then, in the D < case, the radial coordinate is bounded 
as 



/.N r2r-4sn (C, fc) 
r4-r2cn (4, fc) 

/J.-, ri{r3 - r2) + r2{rA - r3)c.n (^, fc) 
7-3 — ^2 -I- {ri — r3)cn (^, fc) 



(ifi) 



/3r4[l-cn(C2, fca)] 
/3 - a - (a -h /3)cn (^2, fc2) 



(A.25) 



r ^ r2, or r^ ^ r. (A. 14) 

If the first inequality holds, we write 

Ta^ri, rt = r2, r^^r^, r^ = r4, (A.15) 

but, if the second inequality holds, we write 



. ^/er3{r2-ri) 
^= 2h '^' 



k = j!jtl^, (A.26) 
r3(r4 — r2) 



where (i) applies if there are 4 real roots and r < 7-2, (ii) is 
used for 4 real roots but r > r^. and (iii) applies when there 
are just two real roots. 



Ta = ra, Th = Ti, re = Vi, r^ = rs. 



(A. 16) 
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